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Abstract 

We extend the model of surface granular flow proposed in jj| to 
account for the effect of an external 'wind', which acts as to dislodge 
particles from the static bed, such that a stationary state of flowing 
grains is reached. We discuss in detail how this mechanism can be 
described in a phenomenological way, and show that a flat bed is lin- 
early unstable against ripple formation in a certain region of parameter 
space. We focus in particular on the (realistic) case where the migra- 
tion velocity of the instability is much smaller than the grains' velocity. 
In this limit, the full dispersion relation can be established. We find 
that the critical wave vector is of the order of the saltation length. Wc 
provide an intuitive interpretation of the instability. 



1 Motivation 

Common observations suggest that flat sand surfaces can become unstable 
when subjected to moving air or water. After some time regular patterns ap- 
pear, as can be observed on desert dunes, underwater sand, 'dry' snow, etc. 
These patterns resemble surface waves; however their physics is completely 
different since in the case of sand there is no surface tension. Following 
Bagnold (H, chap. 11) these patterns can be classified into ripples, ridges 
and dunes. The repetition distance of ridges varies with time, whereas rip- 
ples exhibit a stationary wavelength after some transient. Early qualitative 
arguments by Bagnold ^ suggested that the ripple wavelength A is related 
to the typical path length of the blown grains, called the 'saltation length' 
£. A more quantitative 'two-species' model was proposed by Anderson |J, 
which describes the coupling between the moving grains and the static bed. 
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Such a model predicts that a flat surface is unstable for all wavelengths, 
with a faster growing mode indeed comparable to the typical jump length 
of the grains. However, this model is incomplete: while the dynamics of the 
static bed is treated exactly, the description of the moving phase is highly 
simplified. Alternatively, there are also several numerical models for ripple 
formation Q. In this paper, we extend Anderson's theoretical model of 
ripple formation, by adapting the phenomenological equations for surface 
flow introduced in in the context of avalanches, and further discussed in 

It is worth recalling, after Bagnold ]2| , some basic facts about the motion 
of the grains and the formation of these patterns: (i) There are two quali- 
tatively distinct transport mechanisms for the grains, saltation and surface 
creep Q. The trajectories of grains in saltation is determined by the velocity 
profile of the wind, the air friction limiting the grain velocity and by the 
initial energy of the grain when first expelled from the sand bed. One of the 
characteristic features of the path is the flat angle of incidence which varies 
between 10° and 15°. (ii) The saltation has two effects on the surface: it 
either rebounds and/or ejects grains leading to a new saltation or it pro- 
duces surface creep. There is however no sharp boundary between these two 
processes, since the energy of the ejected grains varies continuously. Both 
saltation and creep lead to a net flow of grains in the direction of the wind, 
(iii) The time scale of ripple formation is much larger than that of salta- 
tion, (iv) The migration velocity of the ripples is much smaller than a mean 
transport velocity (averaged over saltation and creep). 

2 A 'two-species' model with wind 

The phenomenological approach we consider in the following is based on 
the observation that two different species of grains enter the problem: mov- 
ing grains and grains at rest. We will not distinguish between grains in 
saltation and creep, but introduce an appropriately averaged quantity de- 
scribing grains that are convected by either of the two mechanisms, which 
we call the moving grain density R(x,t) 0, where x is the coordinate in the 
direction of the wind and t the time. (We will assume that the problem is 

lr To which one should also add 'suspension', corresponding to very small grains flying 
high in the air. 

2 In principle, one should consider a density R(x, v, t) which depends on the instanta- 
neous velocity of the grains. R(x, i) is the average of R(x, v, t) over all velocities. 
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translationally invariant in the direction transverse to the wind; see || for 
an extension to two dimensions). The grains at rest contribute to the local 
height h(x, t) of the static bed. The dynamical equations for R and h read, 
in the hydrodynamical (long wavelength) limit: 

d t R = -Vd x R + D x e%R + T[R,h\ 

dth = -r[R,h) (i) 

where V and D\ are the average velocity of the grains and the dispersion 
constant, related to the fact that grains do not all move with the same 
velocity []. T describes the rate with which a grain at rest is converted 
into a moving grain (or vice versa) and depends both on R and on the local 
surface profile. For simplicity we have defined R to have the same dimension 
as h, and it can be thought of as the width of grains which has been removed 
from the static bed. Correspondingly, T has the dimension of a velocity. The 
construction of T is based on phenomenological arguments p], and encodes 
different physical processes: 

• Due to the presence of wind, grains can be 'spontaneously' ejected 
from the surface, even in the absence of already moving grains. The rate 
at which this occurs depends on the local wind velocity (or rather velocity 
gradient at the surface). Since the wind velocity tends to be larger when 
the local slope is facing the wind, we write: 

T sp = a + a\d x h - a 2 d x h (2) 

where the coefficients a are positive or zero. We have also included the 
second derivative contribution with a minus sign, since grains are harder to 
dislodge in troughs than at the top of a crest. Note that all these coefficients 
are expected to depend on the external wind velocity. In particular, as shown 
by Bagnold himself, the coefficient «o is only non-zero above a certain critical 
wind velocity, noted Vg uid . 

• When hitting the ground, a moving grain can either be captured or 
transfer a part of its kinetic energy to other static grains and provide new 
moving particles. The rate at which both these process occur is proportional 
to R (at least for small enough R - see below), and also depends on the wind 
velocity and on the local slope. For example, flying grains have a larger 
probability to land on a surface facing the wind, rather than in the wind 

3 These terms can be understood, more generally, as the long-wavelength limit of a 
more general non-local convection term of the kind J K(x — x')R(x' , t) dx' . 
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shadow. This suggests to write the stimulated conversion rate as: 

T st = -R[*y + 7i cU + 72^/*] (3) 

The sign of 70 depends on the strength of the wind; for small wind ve- 
locity, one expects capture to be more important than emission, and thus 
that 70 > 0. As again shown by Bagnold, a localized source of moving 
grains tends to die away when the wind velocity is less than a certain 
impact < ^fl*uid' whereas a steady saltation is found for larger velocities, 
suggesting that 70 < for V > VJ* act . In this case, however, it is easy to 
see that R increases exponentially, and that higher order terms are needed 
to describe the stationary situation. One can think of several non-linear 
effects: for example, collision between flying grains leads to dissipation and 
hence to a poorer efficiency of the impacts on the static bed. Also, the pres- 
ence of a layer of moving grains screens the hydrodynamical flow, which in 
turn reduces the energy transfer between the wind and the saltating grains. 
These effects can be described by a term proportionnal to —f3R 2 in T s t []. 
If trapping dominates (as is the case for under water ripples) one expects 
71 > because more grains fall on the slope facing the convective flow. For 
the same reason, if stimulated emission dominates, as is the case for wind 
blown sand, one expects that 71 < 0. Finally, 72 is positive since, again, 
grains are easier to dislodge at the top of a bump. 

The total conversion rate T is obtained as the sum of T sp and T s t, while 
the model proposed in |lj did not contain the wind induced contribution 
proportionnal to a, nor the non-linear term. The equation for h thus reads: 

d t h = (#70 - a ) + f3R 2 + [Bqx - ai )d x h + {Rj 2 + a 2 )d 2 h (4) 

The gradient term can be interpreted as a translation of the surface profile 
with time, at velocity a% — Rji. The direct action of the wind («i) is indeed 
to erode grains from the windward slope of a bump and transport them in 
the direction of the wind. The other contribution (-R71), however, moves the 
bumps 'backwards' since grains are effectively deposited on the windward 
slope, contributing to a translation of the bump against the wind. (A similar 
discussion can be found in pi [7|, M.) 

Note that the above set of equations is non-linear, so that non-trivial 
dynamics is expected. Some essential features of the model can be inves- 
tigated by linearizing the system in the vicinity of the situation where the 

4 In principle, the dependence of V on R should also be taken into account. We do not 
consider this here, since this does not affect the linear instability analysis. 
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Figure 1: Stability diagram, showing dR/dt as a function of R in an homogeneous 
situation. The case 70 < corresponds to blown sand with V > Vjj^ pact , where 
stimulated emission is very efficient, and where Rq can be non-zero even if oeo = 
(i.e. when V < V^ uid ). The situation where capture dominates (70 > 0) is probably 
relevant for sand under water. 



surface is flat (ho = 0). The moving grain density is then equal to (see Fig 



3 Stability analysis 

We will perform a stability analysis, i.e. investigate whether a small pertur- 
bation is amplified or dies out with time. Therefore we consider R = Rq + R, 
h = ho + h and neglect second order terms of the kind Rh, R? and h 2 . For 
simplicity of notation we drop the bars; the linearized equations then read 



1): 





d t R = -i Q R-Vd x R + D 1 E%R + Wd x h-D 2 6%h + ... 

d t h = ^ R-wd x h + D 2 d 2 x h + ... 



(6) 
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with an effective velocity W = a\ — RqJi and an effective diffusion constant 
D 2 = -R072 + «2 > 0. 70 is equal to 70 + 2(3Rq and is thus always positive. 
A Fourier analysis of the linearized equations leads to 



-iuj 



70 — ikV 
7o 



k 2 D-i 



-iuj 



ikW + k 2 D 2 
ikW - k 2 D 2 




(7) 



where the tilde denotes the Fourier transforms. This system has a non- 
trivial solution if the determinant of the above matrix is zero, leading to the 
relation 

J 1 + uj(a + ib) + (c + id) = 0. (8) 

The coefficients read 



a 
b 

c 
d 



(V + W)k 

- 70 + {D l + D 2 )k 2 

VWk 2 - DiD 2 k 4 
-(DiW + D^k 3 ; 



(9) 



they are functions of the wave vector k and of the system's parameters (V, 
W, Di, D 2 , 70). 

Equation @ establishes a dispersion relation to(k) with two branches 
corresponding to the two solutions of the quadratic equation, where u has 
to be considered as a complex variable. (Writing down the corresponding 
equations for the real and the imaginary part of uj leads to quartic equations.) 
In the context of a stability analysis we are interested in the imaginary part 
of oj{k): as long as it is positive e lu3t will decay exponentially, while a negative 
imaginary part does lead to an instability. This imaginary part is given by: 



2 Im(u; ± ) 
V2 



(a 2 - b 2 - 4c) + (a 2 - b 2 - Ac) 2 + (2ab - Adf 



1/2' 



1/2 
(10) 



which is a function of k. A critical wave vector k* can be defined such that 
Im(o;) is exactly zero, which leads to d 2 — abd+b 2 c = 0. Inserting the explicit 
expressions (S), one finds a cubic equation for k* 2 . Whenever this equation 
admits a positive solution, there will be a finite band of wave vectors [0, k*\ 
which are unstable (see Figure 2). 
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Figure 2: Rescaled damping rate as a function of the rescaled wave vector. The 
plot shows data for r\ = 0.1, 70-D1 = V 2 and D 2 /Di = 0.1 



It is instructive to study the asymptotic behaviour of the functions 
Im((j_). One finds: 



with 77 = W/V. The transport velocity V (by convention) and the diffusion 
constants are positive; the main control parameter remaining is the relative 
migration velocity rj. One sees that for < 77 < 1 there is indeed a band 
of instable wave vectors (one can see from the asymptotic solutions that 
the sign changes for large fc's). The situation where 77 < is stable and 
rj > 1 (i.e. a bump moving faster than the transport velocity) does not seem 
physical. The second branch Im(u;+) is always positive and is thus of no 
importance for our stability considerations. 

Following the intuition that the ripples move much more slowly than the 
grains are transported, we will assume in the sequel that < 77 <C 1, which 
we attribute to the fact that the a coefficients are small compared to V. 




for k <C 7o/F 
for k > 70/^ 



(11) 
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Since D 2 oc a, this suggests that the diffusion constants D\ and D2 are in 
the same ratio, so we write: D2 = SrjDi, where 5 is of the order of one. 
These assumptions make it possible to simplify the algebra and to find the 
solution: 

Im(W - } = 11 7o 2 + fc 2 (2 W + y 2 ) + ^ + ° iV } - (12) 

which is plotted in Figure 2. As we discuss now, three relevant facts can be 
verified with this formula: (i) the critical wave vector is of the order of the 
inverse mean saltation length, (ii) the ripple velocity is of the order of rjV 
and (iii) the time scale of ripple formation is much larger than the saltation 
time scale. 

Let us first give some arguments for (i). Since the saltation trajectories 
result from some random initial vertical velocity of the grains, the saltation 
lengths will also be random, with both short jumps (actually corresponding 
to creep) and long jumps. It is reasonable to assume that the width of the 
saltation length distribution is of the same order as its mean £ (a similar 
assumption is discussed in S). In this situation, the 'Peclet' number defined 
as Pe = V£/D\ is of order one: convective and diffusive effects are of the 
same order of magnitude. In the case where the jump length distribution is 
sharply peaked around £, one would rather have Pe 3> 1. 

Defining £ = Vr, where r is the typical saltation time, one finds that 
the the zero of ( p"2| ) is located at: 

FV ' /(tot + Pe) 2 + 4%t/5 - (Pe + 7o r)| (13) 
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Since 70 is the (renormalized) rate of sticking, it is reasonnable to assume 
that 7oT ~ 1, thereby confirming that k* ~ for Pe ~ 1. On the other 
hand, for weakly dissipative collisions (hard grains) one expects that 70T <C 
1, leading to larger unstable wavelengths ~ £v /( V7o' r - 

The ripple velocity is given by the corresponding dispersion relation, i.e. 
the real part of to{k). One finds: 

The formula shows that for k ~ k*, both phase and group velocities are 
indeed of the order of r\ V, establishing thus (ii). 

Finally knowing the fastest growing wave vector, one finds that the ripple 
formation time t r ippie (determined by the depth of the minimum in figure 2) 
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is a factor I/77 larger than say I/70 or r, i.e. that ripple formation occurs 
on much slower time scales than any microscopic process. The ratio of 
formation time and microscopic time scales should indeed be roughly the 
same as that between migration and convection velocity (hi). 

4 Physical discussion and open questions 

Let us finally give an intuitive interpretation of the instability. Imagine a 
flat surface with a finite number of moving grains above it (i.e. the sta- 
tionary solution). Now imagine a small perturbation of this situation, say a 
small hump. The term dtR ~ d x h in the linearized equations (g) increases 
locally the concentration of the moving grains thus producing a 'cloud' at 
the windward side of the hump. This cloud is convected with the velocity 
Vand after a time unit of I/70 the cloud has moved a distance £ where the 
cloud starts to 'rain' (i.e. moving grains are converted into grains at rest). If 
the position of the hump has in the same time moved in the same direction, 
its height will increase, leading to an instability. (Conversely, if the bump 
moves backward - i.e. if W < - the 'rain' will rather fill the hole and smear 
out the bump). The presence of the diffusive processes counterbalances the 
amplification for small distances and some optimum wavelength of the order 
of £ (corresponding with the minimum in figure 2) becomes visible. 

Summarizing, we have thus shown that equations (1), which are phe- 
nomenological, but motivated by clear physical processes, indeed show an 
instability which is consistent with some essential features of ripple for- 
mation. It is worth noting that our analysis, which concentrated on the 
linearized system in the vicinity of the stationary solution, is universal in 
the sense that a whole class of models behaves in an analogous way (with 
some possible redefinition of the coefficients). For example, a non-linear 
dependence of the velocity V on R does not modify the above analysis, up 
to a redefinition of V. Note also that all phenomenological coefficients are, 
at least in principle, measurable in situations independent from ripple for- 
mation (since they are diffusion constants, convection velocities, deposition 
rates etc.). In this sense it should be possible to check experimentally for 
the consistency of the above description. 

Our conclusions are very similar to those reached by Anderson ||, on the 
basis of a simplified model where the flowing phase (what we have called 
R above) is assumed to be in equilibrium from the outset, and where a 
rather arbitrary distinction is made between 'saltating grains' which are 
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never captured by the bed, and 'reptating' grains which are captured after 
exactly one jump. Correspondingly, the structure of the dispersion relations 
differ in the two approaches. Furthermore, it is difficult to extend Anderson's 
model beyond the linear instability analysis while our model, in principle, 
can account for non- linear effects fl(i| ]. 

Finally, there are several open questions which we would like to mention 
and leave for future work: (i) Can one establish some precise relations be- 
tween the 'microscopic' coefficients (like wind velocity, polydispersity, elas- 
ticity etc.) and the phenomenological parameters? (ii) How are the above 
results modified if one considers two spatial dimensions? Is there an insta- 
bility corresponding to the transverse wavelike shape of the ripples known 
from field observation? (iii) What is the ripple shape and height predicted 
from a non-linear analysis of the equations? (iv) Is it important to consider 
a non-local convection term, rather than the hydro dynamical form written 
in (1)? The question arises since the relevant wavelength is precisely of the 
same order as (and not much larger than) the jump length ^. 
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